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The electronic structure of a vortex line trapped by an insulating columnar defect in a type-II su- 
perconductor is analysed within the Bogolubov-de Gennes theory. For quasiparticle trajectories with 
small impact parameters defined with respect to the vortex axis the normal reflection of electrons 
and holes at the defect surface results in the formation of an additional subgap spectral branch. The 
increase in the impact parameter at this branch is accompanied by the decrease of the excitation 
energy. When the impact parameter exceeds the radius of the defect this branch transforms into the 
Caroli-de Gennes-Matricon one. As a result, the minigap in the quasiparticle spectrum increases 
with the increase in the defect radius. The scenario of the spectrum transformation is generalized 
for the case of arbitrary vorticity. 

PACS numbers: 74.45.+C, 74.78.Na, 74.78.-w 



I. INTRODUCTION 

The study of magnetic and transport properties of 
type-II superconductors in the presence of artificial pin- 
ning centers is known to be an important direction in the 
physics of vortex matter—. The artificial pinning provides 
a unique possibility to control the critical parameters of 
superconducting materials which are important in vari- 
ous applications. For instance, the critical current j c and 
the irreversibility field Hi rr can be enhanced by the in- 
clusion of normal particles and nanorods^, by introduc- 
ing arrays of submicrometer holes^£, and by proton^ and 
heavy-ion irradiation^. Pinning of flux lines appears to 
be especially strong for the case of columnar defects elon- 
gated nearly parallel to the applied magnetic field, when 
vortices can be pinned over their entire length. These 
columnar defects are now widely used to trap vortices 
and to increase the current carrying capacity of super- 
conductors. 

Within the London approximation the interaction be- 
tween a single vortex and an insulating cylindrical cav- 
ity of radius R <C A, where A is the London penetra- 
tion depth, was considered in the pioneering paper— for 
bulk type-II superconductors. For a multiquantum vor- 
tex it was shown 8 - that the maximum number of flux 
quanta which can be trapped by the cylindrical cavity 
is restricted by the value -R/2£, where £ is the super- 
conducting coherence length. The generalization of the 
results of Ref£ for cylindrical cavities of radii R > A has 
been obtained in Ref£. An efficient image method ap- 
propriate for the analysis of the vortex-defect interaction 
in the limit of rather large A values has been developed 
in 10 - 11 . The formation of superconducting nuclei with 
nonzero vorticities near the columnar defects or in per- 
forated films has been studied in Refs i 12 ' 13 . 

Certainly the phenomenological approaches used in 
most of the works cited above can not describe the elec- 
tronic structure of the vortex states in the presence of 
small cavities or colomnar defects of the radius smaller 
than the coherence length £. This issue is closely re- 
lated to the problem of microscopic nature of pinning 



addressed previously in Ref.— for a particular case of 
point-like defects with the scattering cross section much 
smaller than the £ 2 value. An appropriate modification 
of the quasiparticle spectra caused by a single impurity 
atom placed in a vortex core has been studied in^. The 
case of vortices trapped by normal metal cylindrical de- 
fects has been addressed i n 16 i 17 . The interest to micro- 
scopic calculations of electronic structure of the vortex 
states is stimulated by low-temperature scanning tun- 
neling microscopy (STM) experiments which provide de- 
tailed spatially resolved excitation spectr a 18 ' 19 : 20 . The 
modern STM techniques could provide us the informa- 
tion about the number and configuration of the spectral 
branches crossing the Fermi level. Recent STM exper- 
iments on NbSe2 single crystals with a regular array of 
submicron Au antidots have provided images of both sin- 
gle quantum Abrikosov vortices and multiquanta vortex 
states forming near normal antidots^i. 

The goal of our paper is to analyse the transforma- 
tion of the quasiparticle excitation spectra which occurs 
in a vortex pinned by a columnar defect of finite radius 
R % £• We focus on the modification of the anomalous 
energy branches caused by normal reflection of quasipar- 
ticles at the columnar defect boundary. To elucidate the 
key points of the present work we start from the qualita- 
tive discussion of the spectrum transformation scenario. 
Let us consider a vortex pinned at an isolating cylinder 
of a radius R (see Fig. I). The spectrum of quasiparti- 
cle states can be analysed considering one-dimensional 
quantum mechanics of electrons and holes along a set 
of linear quasiclassical trajectories. Each trajectory is 
defined by the impact parameter b and the trajectory 
orientation angle (see Fig.l). For small impact param- 
eters b < R the trajectories experience a normal reflec- 
tion from the defect surface. Hereafter we assume this 
reflection to be specular. Far from the reflection point 
O the superconducting gap is homogeneous (A = Ao) 
and the corresponding superconducting phase difference 
Sip between the trajectory ends is defined by the impact 
parameter b: 5ip = 2 arcsin(|b|/i?). Neglecting the de- 
tails of the inhomogeneous profile of the order param- 
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FIG. 1: Specular reflection of a quasiclassical trajectory at 
the defect surface. 

eter inside the vortex core we can take the gap func- 
tion in the form: A(s) = A exp(i arcsin(|6|/i?) signs), 
where s is the coordinate changing along the trajectory. 
The one-dimensional quantum mechanical problem with 
such order parameter is equivalent to the one describ- 
ing a single mode Josephson constriction^ 2 -. The subgap 
spectrum in this case is known to consist of two energy 
branches: ef(6) = ±A cos(8tp/2) = ±A ay /l - b 2 /R 2 , 
which correspond to the opposite momenta of quasipar- 
ticles propagating along the trajectory. Thus, for small 
impact parameters the scattering of quasiparticles at the 
defect surface is expected to result in the formation of 
new energy branches which arc splittcd from the con- 
tinuum. Taking the quasiclassical trajectories with large 
impact parameters b > R one can see that these trajecto- 
ries are not perturbed by the scattering at the defect and, 
as a consequence, the spectrum in this case should be 
described by the well-known Caroli-de Gennes-Matricon 
(CdGM) expression^. The crossover between two differ- 
ent regimes occurs in the region b ~ ±_R, which should 
be certainly treated more accurately (see below). One 
can expect that in this region the CdGM energy branch 
eo(b) is cut at the energies e ~ ±Eq(R) and transforms 
into the spectral branches £j approaching ±Ao with the 
further decrease in the \b\ value. The resulting spectrum 
as a function of a continuous parameter b does not cross 
the Fermi level: there appears a minigap ~ Eo(R). For 
R <C £ this minigap can be approximately written as: 
6q(R) ~ Ao-R/£. The increase in the defect radius is ac- 
companied by the minigap increase and for R 3> £ all the 
subgap states appear to be only weakly splittcd from the 
±A value. 

The paper is organized as follows. In Sec. [IT] we briefly 
discuss the basic equations used for the spectrum calcu- 
lation. In Sec. [HI] we study the quasiparticle spectrum 
transformation for a singly quantized vortex pinned at a 
columnar defect. In Sec. lIVI we generalize our analysis for 
the case of a multiquantum vortex trapped by the defect. 
We summarize our results in Sec. [Vj 

II. BASIC EQUATIONS 



field B is assumed to be parallel to the cylinder axis 
z. We assume the system to be homogeneous along the 
z— axis, thus, the k z — projection of the momentum is con- 
served. The quantum mechanics of quasiparticle excita- 
tions in a superconductor is governed by the two dimen- 
sional BdG equations for particlelike (it) and holelike (v) 
parts of the two-component quasiparticle wave functions 
i$f(r,z) = (it, v) cxp(ik z z): 

K 2 

- - — (V 2 + kj) it + A(r) v = eu (la) 
2m 

— (V 2 + k 2 ± ) v + A*(r)u = ev . (lb) 
2m 

Here V = d x iiQ + d y yo, r = (x, y) is a radius vector in 
the plane perpendicular to the cylinder axis, A(r) is the 
gap function, and k\ = k 2 F — k 2 . 

Following the procedure described i n 24 i 25 ' 26 we intro- 
duce the momentum representation: 

^=0 = (w/ d2pe,pr/fi|(p) (2) 

where p = |p| (cos9 p ,sin6> p ) = ppo- The unit vector po 
parametrized by the angle 9 p defines the trajectory direc- 
tion in the (x, y) plane. We assume that our solutions 
correspond to the momentum absolute values p close to 
the value hk±: p = hk± + q (\q\ <C hk±). As a next step, 
we introduce a Fourier transformation: 

<MP) = ^ / dse^-WM'faep). (3) 

— oo 

Finally, the wave function in the real space r = 
r(cos 8, sin 0) is expressed from Eqs.([21 in the following 
way (see Ref<2£): 

2tt 

e ^rcos(9 p -9)^ rcos{9p _ e) ,Q p) _JL _ (4) 



where (r, 9, z) is a cylindrical coordinate system. The 
boundary condition at the surface of the insulating cylin- 
der requires 

^R,9) = ( U ) = (5) 

\V / r=R 

— / d9 p c lk±Rcos( - e ?- d) ip(Rcos(9 p - 9),9 P ) = 0. 

27T J 


To obtain the Andrccv equations along the trajectories 
we look for a solution in the eikonal approximation 

$(s,6 p ) = e iS ^g(s,6 p ) 

assuming g to be a slowly varying function of 9 P . Quasi- 
particles propagating along the classical trajectories par- 
allel to = k±(cos9 p ,sm8 p ) are characterized by the 
angular momentum fi = —k±b, where 



Hereafter we consider a columnar defect as an insu- 
lating infinite cylinder of the radius R. The magnetic 



(0) 
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is the trajectory impact parameter. Assuming the vor- 
tex axis to coincide with the cylinder axis we obtain the 
axially-symmctric problem with the conserved angular 
momentum fi. 

Finally, the quasiclassical equations for the envelope 
g(s, Op) read: 



dg 

ifiVx& z — + & x ReA(r)g 

OS 



(7) 



o-j / ImA(r) g = [ e + \x 



where &i are the Pauli matrices, mV± = hk±, loh = 
\e\H/mc is the cyclotron frequency, and 

x = s cos 6 p — b sin 9 p , y = s sin 9 p + b cos 9 p , 
x ± iy = (s ± ib) e ±ie * . 

The term proportional to luh can be included to the en- 
ergy as an additive constant (see also21): 

huj H 

e = e+—». 

Our further analysis of quasiparticlc excitations is based 
on the Andrecv equations (JT]) which must be supple- 
mented by the boundary condition 



III. SINGLY QUANTIZED VORTEX PINNED 
BY A COLUMNAR DEFECT 

We now proceed with the analysis of the subgap spec- 
trum for a singly quantized vortex trapped by the colum- 
nar defect of the radius R. The order parameter A(x,y) 
takes the form 



A = A S v {r)c lf> , r=y/x 2 + y 2 >R. 



(8) 



Here 5 V (r) is a normalized order parameter magnitude for 
a vortex centered at r = 0, such that S v (r) = 1 for r — > 
oo. In (s, Op) variables one obtains for r = Vs 2 + b 2 > R: 



A = D b (s) e ie * , D b (s) = A $vW f \P (s+*b) . (9) 

The cylindrical symmetry of our system allows to sepa- 
rate the Op— dependence of the function g: 



g( Sl O p ) = e* a ^ 2 f(s) 



(10) 



The total wave function if)(s,O p ) should be single valued 
and, thus, the angular momentum fi is half an odd inte- 
ger. The quasiclassical equations ([7]) take the form 



ihV ± a z dJ + A b (s)f = ef. 



(11) 



is the gap operator. Changing the sign of the coordi- 
nate s one can observe a useful symmetry property of 
the solution of Eq. (flT|) : 

f(-s) = Ca y f(s) , (13) 
where C is an arbitrary constant. 

A. Boundary condition. 

As a next step we rewrite the boundary condition ([5|) 
for wave functions /(s) defined at the trajectories. Re- 
placing Op by a — Op — and shifting the limits of inte- 
gration in Eq.((5|) we find: 



2n 



da e 



;k±R cos ct+ipOL 



* /2 f(R cos a) 



= 0. 



(14) 



Assuming kj_R 3> 1 and the function e l<Jza / 2 f(s) to vary 
slowly at the atomic length scale we evaluate the above 
integral using the stationary phase method. For a given 
value of angular momentum the stationary phase points 
are given by the condition: sin 0^2 = n/k±R = —b/R. 
One can see that for |6| > R the stationary phase points 
disappear and, result, the integral p4[) is always 

vanishingly small. In this case the boundary condition 
at the cylinder surface does not impose any restrictions 
on the wave function / defined at the trajectories. In the 
opposite limit \b\ < R one can find two stationary angles 
Q'i = ao = — arcsin(6/i?) and 02 = tt — ao which are in 
fact the orientation angles for an incident and specularly 
reflected trajectories shown in Fig. [T] Summing over two 
contributions we can rewrite the boundary condition (fTJ 
as follows: 



e i * 1 f(s ) = e- i ^f(-s ), 



(15) 



where 



where Sq = VR 2 — b 2 , 2/3o = ao — tt/2 and 
fix = kj_s + (2/j, + <T 2 )/?o - 37r/4 . 

B. Solution for large impact parameters |6| > R. 

In this case the quasiparticle states at the trajecto- 
ries are not affected by the the normal scattering at the 
columnar defect boundary and the behavior of an anoma- 
lous energy branch is described by the standard CdGM 
solution for a single Abrikosov vortex. For the sake of 
completeness we give below the expressions for this spec- 
trum and the corresponding wave functions. 

Let us follow the derivation in Ref3i and consider the 
imaginary part of the gap operator (|12|) as a perturba- 
tion. Neglecting this term in Eq. (jlip we find: 



A b (s) = a x ReD b (s) - a y lmD b (s) 



(12) 



ihV±a z d s f + a x ReD b (s) f = sf 



(16) 
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The above equation has a zero eigenvalue e = with 
the following expression for the corresponding normalized 
eigenfunction / : 



(17) 



where 



S +OC 

K (s) = j^J dtReD b (t) , 7 = / d S e~ 2K °^ . 

-co 

(18) 

The first order perturbation theory gives us the CdGM 
excitation spectrum So(b) for |&| > R: 



1q J Vs + b z 



which occurs at b ~ ±R. To develop an analytical de- 
scription of this crossover we choose to apply the method 
used above to derive standard CdGM expressions and 
based on the perturbation theory with respect to the 
imaginary part of the gap function. One can expect this 
method to be most adequate for the crossover region of 
b ~ R. As for the limit b < R we shall check the valid- 
ity of this method using the comparison with our direct 
numerical analysis of Eq. (|2"Tj) . 

Neglecting the maginary part of G we find an exact 
solution of the equation (f2"Tj) corresponding to zero energy 
£ = 0: 



K(s) 



(24) 



where 



s +oo 

K{s) = j^p- J dtRcG(t), 1= Jdse- 2K{s) , (25) 



C. Solution for small impact parameters |6| < R. 

In this case the specular reflection at the cylinder sur- 
face changes the trajectory direction and strongly modi- 
fies the spectrum. The boundary condition at the surface 
r = R is determined by the equation (| 1 5[) . Let us intro- 
duce the function 



F(s) 



e +i ^/(s + so), s>0 
e-^fis-so), s<0 



(20) 



which is defined at the full s axis and appears to be con- 
tinuous at s = 0: F(-0) = F(+0). The equation for F 
reads: 



ihV ± a z d s F + fr x Re G(s) F 

- fry Jm G(s) F = eF 



(21) 



where 



-A 



W(N + so) 2 + fr 2 ) 
V(N + so) 2 + 6 2 



sb/R + i[ R + \s\y/l -b 2 /R 2 



G(s) 



Taking the limit of large |s| we find: 

G(s) = -iA e iQo|s|/s 



(22) 



(23) 



One can see that in agreement with the qualitative ar- 
guments given in the Introduction the phase difference 
between the opposite ends of the trajectory equals to 
Sip = — 2«o = 2 arcsin(6/i?). Provided we neglect the 
inhomogeneity of the order parameter phase inside the 
core we find the resulting expression for the spectrum: 
e = ±Ao \Jl — b 2 / R 2 . Certainly, such simplification does 
not allow us to study the crossover to the CdGM branch 



and x = sign 6. The solutions ([24]) . (f25j) appear to decay 
both at negative and positive s and, thus, we get a local- 
ized wave function describing a bound state. Using this 
localized solution as a zero-order approximation for the 
wave function the spectrum can be found within the first 
order perturbation theory. Note, that our perturbation 
procedure fails for | b \ — » because of the increase in the 
localization radius of the wave function (|24| . 

The first-order approximate solution of the quasiclas- 
sical equations (|2Tj) takes the form: 

F(s) = A f M e -*M +B(s) (J J c K{ - s \ (26) 



where 



B{s) = w: I dt ^ £ - x Im °w e ~ 2K{t) ■ w 

— oo 

To avoid the wave function divergence we should put 
XImG(t)] c~ 2K ^ =0. 



dt 



This condition gives us the excitation spectrum e s as a 
function of the impact parameter b for |6| < R : 



s s (b) = 



XA 



+ 00 



ds 



5vW{\s\+s ) 2 +b 2 ) 
^(\ s \ + S0 ) 2 + b 2 

R + \ s yi-b 2 /R 2 ) C - 2K ^ 



(28) 



It is evident that e s (R) = e (R) and, thus, the expres- 
sions (fT!?]) and (|2^|) describe the spectrum e(b) for an 
arbitrary impact parameter b : 



e(b) = 



£o(b) , 



b\<R 
\b\>R 



(29) 
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The discontinuity of the derivative dejdb at | b | = R ap- 
pears because of the breakdown of the above quasiclas- 
sical description for the rectilinear trajectories touching 
the surface of the defect. 

In figures 2 and 3 we compare the typical plots of quasi- 
particle spectra obtained analytically, i.e., using Eq. ([2^|) . 
and numerically. The qualitative behavior of the spec- 
trum is weakly sensitive to the concrete gap profile inside 
the core and, thus, we choose a simple model profile: 



(30) 



where the core size equals to the coherence length 
£ = S-Vf/Aq. We plot here only the spectrum for pos- 
itive energies and k z momenta because the eigenvalues 
for e < and k z < can be found using the spec- 
trum symmetry properties: e(—b, k z ) = —e(b, k z ) and 
e(b, — fc 2 ) = e(b, k z ). To find the spectral branch e(b, k z ) 
numerically we solve quasiclassical equations pip for 
s > so requiring the decay of the wave function / at 
s — > oo. An appropriate boundary condition for elec- 
tron /„ and hole /„ components of the wave function 
/ = (/«> fv) at s = can be found from Eq. (JTSJ) and the 
symmetry property (|13|) : 



fv(so) 



7u(s ). 



(31) 



For \b\ > R we put So — 0, and the boundary condition 
(f3T|) takes the form /„(0) = -«x/„(0). 

Comparing the spectrum (j2"9"|) with the branches ob- 
tained from the direct numerical analysis of the eigen- 
value problem (fTTj) . (f3"Tj) one can see that the pertur- 
bation method provides a reasonable description of the 



3-0-0-0 




-0.2 0.0 0.2 0.4 0.6 0.8 1.0 



FIG. 2: (Color online) The spectral branches as functions of 
the impact parameter b, obtained from numerical solution of 
the eigenvalue problem (|ll|l . (|31|l . are shown by red triangles 
and circles for k z = and k z — 0.9&F , respectively. The 
spectral branches calculated using Eq. (|29[) are shown by blue 
dash lines. Here we put R — 0.1£. 
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FIG. 3: (Color online) The quasiparticle spectra obtained 
from the numerical solution of the eigenvalue problem 
(|31[) are shown by red circles (n = —10.5), triangles (n — 
— 18.5), and squares (/i = —25.5). The spectral branches 
calculated using Eq. (129[) are shown by blue dash lines. Here 
we put R = 0.1£, fcrC = 200. 



energy spectrum behavior in a wide range of the impact 
parameters. As one would expect, the perturbation pro- 
cedure fails for small impact parameters \ b\ <C i?. Con- 
trary to the CdGM case the spectrum branch (|29|) does 
not cross the Fermi level, and the minigap in the quasi- 
particle spectrum A m j n = s(R) grows with the increase 
in the cylinder radius R (see Fig. 2). Existence of the 
minigap in the spectrum of quasiparticles should result 
in peculiarities of the density of states (DOS) and can be 
probed by the STM measurements. For | /i | < kpR the 
spectrum e{k z ) has a minimum (see Fig. 3), therefore we 
can expect the appearance of a van Hove singularity in 
the energy dependence of the DOS. 



IV. QUASIPARTICLE SPECTRUM OF A 
MULTIQUANTUM VORTEX 

In this section we generalize the above analysis for 
the case of a multiquantum vortex pinned by the colum- 
nar defect of the radius R. The multiquantum vortices 
can be trapped at columnar defects either for a rather 
large defect radius or for the mixed state in mesoscopic 
sample a 21 ' 29 ! 30 ! 31 . In the absence of defects the spectrum 
of a multiquantum vortex with the vorticity M is known 
to consist of M anomalous energy branches 2 ^. The be- 
havior of these branches has been previously investigated 
both numerically and analyticall y 16 ' 17 ' 26 ! 32 ' 33 . Here we 
restrict ourselves by the numerical solution of the eigen- 
value problem (jlip assuming that the order parameter 
Ajf(r) takes the form 



A M (r) = A [S v (r) 



{ M e we 



r> R, 



(32) 
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FIG. 4: (Color online) The spectral branches as functions of the impact parameter b, obtained from the numerical solution of 
the eigenvalue problem ([TTJ, ([36} for M = 2 (a) and M = 3 (b) (R = 0.1£) . 



where the function S v (r) is determined by the expression 
(|30l) . In (s, 6p) variables one obtains for s > sq: 



(33) 



D M (s) = A 



Vs 2 + b 2 



-i M 



(s + ib) M . (34) 



Using the transformation 



(35) 



we can rewrite the quasiclassical equations ([7]) in the form 
(fTTj) with the gap operator 



Aft(s) = a x ReD M (s) - a y lmD M (s) . 



(36) 



The symmetry properties of both the gap operator A&(s) 
and the equation pip depend on the vorticity M: 



A 6 (-s) 



for even M 



-A*(s), for odd M 



This fact allows to obtain the following condition: 
Ca x f(s), for even M 



/(-*) 



C& y f(s), for odd M ' 



(37) 



(38) 



which generalizes the condition (|13p for a multiquantum 
vortex. Here C is an arbitrary constant. Using the sta- 
tionary phase method we can write the boundary condi- 
tion for wave functions f(s) at the surface of the insulat- 
ing cylinder in the form: 



where 



e i ? M f(s ) = e- i ^f(-s ), 



0m = k±s + (2/i + Ma z )(3 - 3tt/4 . 



(39) 



Taking into account the Eq. (|38l) the boundary condition 
(|39|) can be written for electron f u and hole f v compo- 
nents of the wave function /: 



f v (s )=±e lMa °f u (s ) 



(40) 



For | b | > R we can put here so = and «o = — 7r/2. The 
choice of the sign in ([4^]) depends on the number of the 
spectral branch. The typical plots of quasiparticle spec- 
tra obtained from numerical solution of the eigenvalue 
problem (TTTj) . (|36|) with the boundary condition (|40j) for 
vortices with winding numbers M = 2, 3 are shown in 
Fig. 4. Similarly to the case of a singly quantized vortex 
the small b part of the spectrum is formed by the spectral 
branches induced by the normal scattering at the defect. 
These branches transform into the standard anomalous 
ones with the increase in the | b | value. With the increase 
in the cylinder radius all the spectral branches appear to 
be expelled from the Fermi level. 



V. SUMMARY 

To sum up, we described a transformation of the sub- 
gap spectral branches of quasiparticle excitations in vor- 
tices pinned by columnar defects of finite radii. We find 
that the normal scattering at the defect surface results 
in the appearance of additional spectral branches which 
transform into the CdGM one with an increase in the im- 
pact parameter of quasiparticle trajectories. The increase 
in the defect radius is accompanied by the increase in 
the minigap in the spectrum which can be observed, e.g., 
in the STM measurements. One can expect that such 
changes in the spectrum behavior should affect strongly 
the dynamic mobility of vortices in the presence of ac 
transport current (see, e.g.f^l for review). 
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